% This function returns the signal values on a given time resolution. It also
% generates the noise on the same grid.
% the user can specify either the required time resolution (in sec) or the
% required number of points in the signal duration (in which case TimeRes=-1).

function [SignalAmp, noise_unscaled, TimeAxis] = GenerateSignalValues(Signal,MWC,TimeRes,Npts,handles)

% define the time-axis
% we don't plot the signal on its entire duration since it will consume too
% much time...

%globals
global printArray;

TimeAxis = 0:TimeRes:((Npts-1)*TimeRes);  
dt = TimeRes;

% read signal parameters
N       = getfield(Signal.Structure, 'N');
B       = getfield(Signal.Structure, 'B');
fmax    = getfield(Signal.Structure, 'fmax');
SignalType = getfield(Signal.Structure, 'Type');
Carriers = getfield(Signal.Values, 'Carriers');

% generate signal amplitudes
SignalAmp = zeros(size(TimeAxis));

if (strmatch(SignalType, 'sinc'))
        Energy = getfield(Signal.Values, 'SincEnergies');
        Delay = getfield(Signal.Values, 'SincDelays');
        for n=1:(N/2)
            SignalAmp = SignalAmp+sqrt(Energy(n)) * sqrt(B)*sinc(B*(TimeAxis-Delay(n))) .* cos(2*pi*Carriers(n)*(TimeAxis-Delay(n)));
        end
        
elseif (strmatch(SignalType, 'bpsk'))
        Tsymbol = getfield(Signal.Structure, 'BPSK_Tsymbol');
        NumSymbols = floor((TimeAxis(end)-TimeAxis(1))/Tsymbol);
        BitStream = randsrc(N/2,NumSymbols);
        Energy = getfield(Signal.Values, 'BPSKEnergies');
        for n=1:(N/2)
                xBB = zeros(size(TimeAxis));
                for symInd = 1:NumSymbols   % can improve by creating only delat[n] and filtering with the sinc or rcosine functions
                    symPoint = (symInd-1)*Tsymbol;
                    xBB = xBB + BitStream(n,symInd)*sinc((TimeAxis-symPoint)/Tsymbol);
                end
                SignalAmp=SignalAmp+sqrt(2*Energy(n)*B)*xBB.*cos(2*pi*Carriers(n)*TimeAxis);
        end
        
elseif (strmatch(SignalType, 'qpsk'))
        %  We compute the signal symbols only now according to the specified number of points, 
        %  which also dictates the signal duration
        Tsymbol = getfield(Signal.Structure, 'QPSK_Tsymbol');
        NumSymbols = floor((TimeAxis(end)-TimeAxis(1))/Tsymbol);
        Energy = getfield(Signal.Values, 'QPSKEnergies');
        BitStreamI = randsrc(N/2,NumSymbols);
        BitStreamQ = randsrc(N/2,NumSymbols);


        for n=1:(N/2)
                xBBI = zeros(size(TimeAxis));
                xBBQ = zeros(size(TimeAxis));
                for symInd = 1:NumSymbols   % can improve by creating only delat[n] and filtering with the sinc or rcosine functions
                        symPoint = (symInd-1)*Tsymbol;
                        xBBI = xBBI + BitStreamI(n,symInd)*sinc((TimeAxis-symPoint)/Tsymbol);
                        xBBQ = xBBQ + BitStreamQ(n,symInd)*sinc((TimeAxis-symPoint)/Tsymbol);
%                         if mod(symInd,10)==1
%                             updateProgBar((n-1)*NumSymbols+symInd,handles,N*NumSymbols);
%                         end
                end
                SignalAmp=SignalAmp+sqrt(2*Energy(n)*B)*(xBBI.*cos(2*pi*Carriers(n)*TimeAxis)+xBBQ.*sin(2*pi*Carriers(n)*TimeAxis));
        end
elseif (strmatch(SignalType, 'cosine')) %ETGAR
        Energy = getfield(Signal.Values, 'CosineEnergies');
        Delay = getfield(Signal.Values, 'CosineDelays');
        for n=1:(N/2)
            % %plain cosine:
            SignalAmp = SignalAmp+sqrt(Energy(n)) * sqrt(B)*cos(B*(TimeAxis-Delay(n))).* cos(2*pi*Carriers(n)*(TimeAxis-Delay(n)));          
            CurrentCosFreq = B/(2*pi)+Carriers(n);
            save ('Data/CurrentCosFreq.mat','CurrentCosFreq');
            % % cosine*exp(-(Factor)*t) - use with N=2
            %Factor=B/16;
            %SignalAmp = SignalAmp+sqrt(Energy(n)) * sqrt(B)*exp(-(Factor)*(TimeAxis-Delay(n))).*cos(B*(TimeAxis-Delay(n))).* cos(2*pi*Carriers(n)*(TimeAxis-Delay(n)));          
            % % Optional - playing simple music - Doesn't work yet:
            %SignalAmp = SignalAmp+sqrt(Energy(n)) * sqrt(B)*SimpleMusicSig((TimeAxis-Delay(n))).* cos(2*pi*Carriers(n)*(TimeAxis-Delay(n)));    %Doesn't work yet       
        end
        %plotting Cosine:
%         TsCosine=TimeAxis(2)-TimeAxis(1);
%         FsCosine=1/TsCosine;
%         [Cosine_vec,Cosine_freq] = FreqDomain( FsCosine , SignalAmp, fmax);
%         figure;
%         plot(Cosine_vec,(abs(Cosine_freq)));
        
end

%% DEBUG DISPLAY CODE TO EXPLORE THE SIGNAL FOURIER TRANSFORM

if(printArray.signal)
    Fs = 1/dt;                    % Sampling frequency
    T = 1/Fs;                     % Sample time
    L = length(SignalAmp);                     % Length of signal
    t = (0:L-1)*T;                % Time vector
    y=SignalAmp;
    NFFT = 2^nextpow2(L); % Next power of 2 from length of y
    Y = fftshift(fft(y,NFFT)/L);
    f = Fs/2*linspace(-1,1,NFFT);
%     Plot single-sided amplitude spectrum.
    figure,plot(f,2*abs(Y)) 
    title('Single-Sided Amplitude Spectrum of y(t)')
    xlabel('Frequency (Hz)')
    ylabel('|Y(f)|')
    %Adding fp rubrices
    fp=MWC.fp;
    fpVec = [-fp/2:-fp:-Fs/2 , fp/2:fp:Fs/2];
    fpInd = min(Energy)*ones(size(fpVec));
    hold on;
    stem(fpVec,fpInd,'*g');
end


%% Generate noise 
% Load rand generator seed
rng(Signal.Values.SavedRandomStateForNoiseGen);

fnyquist = 2*fmax;
if (fnyquist < 1/dt)  % the time resolution is higher than the Nyquist rate, so we need to reject the noise beyond fmax
    noise_wideband = randn([1,length(TimeAxis)+3500/2] );  %% the extensio is to avooid the delay caused by the filter
    noise_unscaled_temp = filter(fir1(3500,2*fmax/(1/dt)),1,noise_wideband);
    noise_unscaled = noise_unscaled_temp(3500/2+1:end);
else
    noise_unscaled = randn(size(TimeAxis));
end

